arXiv:1509.06894vl [astro-ph.SR] 23 Sep 2015 


MNRAS 000, 000-000 (0000) 


Preprint 24 September 2015 


Compiled using MNRAS style file v3.0 


Validation of solar-cycle changes in low-degree helioseismic 
parameters from the Birmingham Solar-Oscillations 
Network 

R. Howe^*, G.R. Davies^, W.J. Chaplin^ Y.R Elsworth\ and S.J. Hale^ 

^School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom 


24 September 2015 


ABSTRACT 

We present a new and up-to-date analysis of the solar low-degree p-mode parame¬ 
ter shifts from the Birmingham Solar-Oscillations Network (BiSON) over the past 22 
years, up to the end of 2014. We aim to demonstrate that they are not dominated by 
changes in the asymmetry of the resonant peak profiles of the modes and that the pre¬ 
viously published results on the solar-cycle variations of mode parameters are reliable. 
We compare the results obtained using a conventional maximum likelihood estimation 
algorithm and a new one based on the Markov Chain Monte Carlo (MCMC) technique, 
both taking into account mode asymmetry. We assess the reliability of the solar-cycle 
trends seen in the data by applying the same analysis to artificially generated spec¬ 
tra. We find that the two methods are in good agreement. Both methods accurately 
reproduce the input frequency shifts in the artificial data and underestimate the am¬ 
plitude and width changes by a small amount, around 10 per cent. We confirm earlier 
findings that the frequency and line width are positively correlated, and the mode 
amplitude anticorrelated, with the level of solar activity, with the energy supplied to 
the modes remaining essentially unchanged. For the mode asymmetry the correlation 
with activity is marginal, but the MCMC algorithm gives more robust results than 
the MLE.The magnitude of the parameter shifts is consistent with earlier work. There 
is no evidence that the frequency changes we see arise from changes in the asymmetry, 
which would need to be much larger than those observed in order to give the observed 
frequency shift. 


1 INTRODUCTION 

The variation of solar acoustic mode (p-mode) parameters 
with solar activity is one of the best-known results in he¬ 
lioseismology, and has been studied at a wide range of spa¬ 
tial and temporal scales over the last few decades. In the 
light of recent developments in mode parameter estimation 
techniques it seems timely to revisit these findings. In this 
work we re-analyse the Sun-as-a-star data from the Birm¬ 
ingham Solar-Oscillations Network (BiSON) over the last 
two solar cycles and verify our ability to detect subtle vari¬ 
ations using two different algorithms, one that uses conven¬ 
tional maximum-likelihood estimation and a new one based 
on Bayesian principles 


1.1 Historical Background 

1.1.1 Frequency Variation 

Small changes in the frequencies of the low-degree modes, 
positively correlated with proxies for the solar activity such 
as sunspot number and 10.7 cm radio flux, were reported for 
example by Woodard & Noyes (1985) using data from the 
Active Cavity Radiation Monitor (ACRIM), and later con¬ 


firmed in ground-based observations by Palle et al. (1989) 
and Elsworth et al. (1990), while Libbrecht & Woodard 
(1990) first reported solar-cycle changes in the frequencies of 
medium-degree modes from resolved-Sun observations at the 
Big Bear Solar Observatory. Libbrecht & Woodard (1990) 
also pointed out that the variation increased with frequency, 
following the inverse of the so-called ‘mode inertia.’ This in¬ 
dicates that the mechanism responsible for the shifts is lo¬ 
cated close to the solar surface. The frequency dependence 
of the variation is harder to detect in low-degree data due 
to the smaller number of modes available, but within a few 
years it was reported by Elsworth et al. (1994). 

In resolved-Sun data the frequency shifts follow the 
activity variation in location as well as in time (e.g. 
Howe, Komm & Hill 2002, and references therein); modes of 
the same radial order n and degree I but different azimuthal 
order m have different distributions of power with latitude 
and so are differently affected as the magnetic activity belts 
move in latitude over the solar cycle. Even for low-degree 
modes, Chaplin et al. (2004); Jimenez-Reyes et al. (2004); 
Chaplin et al. (2007); Salabert et al. (2015) have reported 
marginally significant differences in the response of the fre¬ 
quency to activity for modes of different degree, consistent 
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with the different latitudinal distribution of power for the 
different spherical harmonics. 


1.1.2 Power and width variations 

Magnetic activity influences not only the frequency of 
the acoustic modes but also their power and lifetime. 
Elsworth et al. (1993) reported an anticorrelation between 
solar activity and the amplitude of low-degree modes, with 
about a 30 per cent change between solar minimum and so¬ 
lar maximum. The positive correlation between mode width 
and activity took longer to establish, due to the difficul¬ 
ties introduced by systematic effects such as the broaden¬ 
ing of the modes when the duty cycle is low. However, 
Chaplin et al. (2000), after a careful analysis of BiSON 
data, reported a change of around 24 per cent in the line 
width and a similar decrease in the mode power (amplitude 
squared) between solar minimum and maximum. The com¬ 
bination of these two results suggested that the solar-cycle 
change affects the mode damping rather than the excita¬ 
tion, with the rate of energy supplied to the modes remain¬ 
ing constant. Komm, Howe & Hill (2000a,b) found activity- 
correlated changes in the mode width for medium-degree 
Global Oscillation Network Group (GONG) data, which 
could also be localized to the latitudes where magnetic activ¬ 
ity is present. (Komm, Howe & Hill 2002). Again, the energy 
supply rate appears to remain constant. 

Salabert et al. (2007) studied the line width variation 
in 9.5 years of BiSON data and found marginal evidence for 
a degree-dependence of the sensitivity, with the I = 0 shifts 
about half the size of those for 1^1, consistent with what 
would be expected from the latitudinal distribution of the 
modes in relation to the latitudinal migration of the activity 
belts during the solar cycle. They concluded that previous 
estimates of Sun-as-a-star line width changes may have been 
overestimated by about 50 per cent because the shifts were 
averaged over all the modes. 


1.1.3 Asymmetry 

As discussed below in Section 2, the mode peaks show 
a small asymmetry due to the correlation of mode exci¬ 
tation with the noise that drives it (Duvall et al. 1993). 
This term is negative for velocity observations and gener¬ 
ally positive for measurements made in intensity; it is small¬ 
est where the modes are strongest and for low-degree ob¬ 
servations can reach a value of a few per cent at the ex¬ 
tremes of the p-mode band and a fraction of one per cent 
at 3 mHz. Jimenez-Reyes et al. (2007) reported a fractional 
change of about 15 per cent between solar maximum and 
solar minimum in the asymmetry of low-degree modes in 
data from BiSON and the Global Oscillations at Low Fre¬ 
quencies (GOLF) instrument on the Solar and Heliospheric 
Observatory (SOHO) spacecraft, with the strongest (most 
negative) asymmetry at solar maximum. As the uncertainty 
in the asymmetry depends strongly on the signal-to-noise 
ratio of the modes and also on the duty cycle, the result was 
most clearly seen in the GOLF data and was described as 
‘marginally significant’ for BiSON. 


1.2 Mechanisms 

There is still some uncertainty about the precise mech¬ 
anisms responsible for the parameter shifts. Attempts 
have been made to model the frequency shifts by invok¬ 
ing the direct effects of magnetic field at the tachocline 
(Roberts & Gampbell 1986), the sunspot anchoring zone 
around 50 Mm below the surface (Foullon & Roberts 2005), 
the photosphere (Bogdan & Zweibel 1985), or the chromo¬ 
sphere (Gampbell & Roberts 1989; Jain & Roberts 1994), 
but none of these have predicted shifts of the observed mag¬ 
nitude. The shifts have also been considered as an indirect 
effect of temperature changes associated with the activity 
belts (Kuhn 1988) and as an effect of a change in acous¬ 
tic cavity size (Dziembowski & Goode 2005). None of these 
models is universally accepted. In any case, it is possible to 
consider the frequency shifts as a solar-cycle proxy in them¬ 
selves, and one of the few that has some sensitivity to layers 
below the photosphere. 


1.3 Motivation 

Two relatively recent developments have prompted us to 
revisit the mode parameter variations over the last two 
decades of BiSON observations. 

First, although there is broad consensus on the size and 
sign of the different solar-cycle effects on mode parameters, 
Korzennik (2013) has reported results from a re-analysis of 
GONG, Michelson Doppler Imager (MDI), and Helioseismic 
and Magnetic Imager (HMI) medium-degree data showing 
frequency changes smaller by about a factor of two than have 
been seen in the standard analysis of these data sets. The 
author attributes the discrepancy to the use of asymmetric 
peak profiles in the fitting, and also reports variations in the 
peak width and asymmetry, with the asymmetry changes 
possibly accounting for the ‘missing’ frequency shifts. These 
results make it timely to re-investigate the frequency varia¬ 
tions over the past two solar cycles of BiSON observations 
and assess the validity of our mode-parameter variations us¬ 
ing fits that take into account the asymmetry of the modes. 

Secondly, there is growing evidence of, and interest in, 
temporal variations in the mode frequency that cannot be 
described simply as a linear function of the global activity 
level. These may take the form of fluctuations superimposed 
on the linear trend with activity proxy or of changes in the 
strength of the trend over time; the two are not necessarily 
clearly distinguished, as the presence of an additional short¬ 
term variation may reduce the strength of the correlation 
with activity. Following the long, anomalous solar minimum 
after Solar Gycle 23 and the relatively weak Gycle 24, there 
has been particular interest in any changes in frequency 
shifts that might indicate changes in the structure of the 
outer solar layers. Basu et al. (2012) found that the shifts in 
Gycle 23 BiSON data did not follow the trend with activ¬ 
ity extrapolated from the Gycle 22 data, with the mismatch 
showing up first in the low-frequency shifts and only later 
at higher frequencies. On the other hand, Salabert et al. 
(2015), using GOLF data, reported a frequency increase in 
Gycle 24, particularly for I = 1 modes, that was larger than 
that expected from the increase in activity level when com¬ 
pared to the results from Gycle 23. It seems that the inter- 
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pretation of such results is quite sensitive to which modes 
are averaged and what period is taken as the reference. 

There have also been reports of variations in the mode 
frequencies on time scales that do not correspond to the 11- 
year solar cycle. Howe et al. (2006) reported aperiodic fluc¬ 
tuations in individual low-degree mode freqencies after the 
subtraction of the activity effect, which were correlated be¬ 
tween different instruments but not with any activity index 
and not correlated between different modes. These varia¬ 
tions were attributed to stochastic fluctuations in the solar 
spectrum. Fletcher et al. (2010); Broomhall et al. (2011); 
Simoniello et al. (2012, 2013) later reported a quasi-biennial 
variation in the mode frequencies, superimposed on the 11- 
year variation, which they interpreted as evidence of a ‘sec¬ 
ond solar dynamo’ located in the near-surface layers. This 
variation appears to be correlated with other activity prox¬ 
ies in high-activity periods but persists during solar minima 
when the proxies are flat. Along similar lines, Jain et al. 
(2011) found that the GONG medium-degree frequencies 
showed hints of a ‘double minimum’ between Solar Cycles 23 
and 24, with the frequencies first reaching their lowest value 
ahead of the activity minimum; this may reflect a similar 
phenomenon to that seen in the low-degree data. 

All of these results make it timely to re-examine the 
parameter variations and assess their reliability. 

1.4 Arrangement 

The arrangement of the paper is as follows: in Section 2 
we discuss the issues related to correlated and uncorrelated 
noise and peak asymmetry and go on to describe the specifics 
of the data simulation. In Section 3 we describe the algo¬ 
rithms used for extracting the mode parameters. In Section 4 
we describe the tests carried out on the artificial data and 
discuss the results and their implications for the analysis of 
BiSON data. In Section 5 we present the results from fitting 
the BiSON data, and in Section 6 we discuss our conclusions. 


2 GENERATION OF ARTIFICIAL DATA 

In order to test for possible biases introduced into the mode 
parameters and their variations we have carried out a se¬ 
ries of tests on artificial data. The generation of the arti¬ 
ficial data was based on the solarFLAG simulator, which 
was used to make data for the hare-and-hounds study on 
rotational frequency splittings discussed by Chaplin et al., 
(2006) and also to test peak-bagging of Sun-as-a-Star he¬ 
lioseismic data in Jimenez-Reyes et al. (2007). The original 
solarFLAG simulator did not include effects of correlated 
excitation, or of correlations of the excitation with back¬ 
ground noise. An important consequence was that the arti¬ 
ficial mode peaks in the frequency power spectrum showed 
no asymmetry, unlike their real solar counterparts. Since 
any analysis which seeks to extract accurate estimates of 
the frequencies must cope with this asymmetry it was felt 
we needed a simulator that could provide such a test. The 
departure of the mode shape from the pure Lorenztian 
predicted by a simple harmonic oscillator has three main 
causes. They are the localization of the noise source, the im¬ 
pact of correlated noise on mode excitation, and finally the 
more subtle effect of the presence of the wings of nearby 


correlated modes. We have used a simple but very pow¬ 
erful method to introduce in the time domain the effects 
of asymmetry, which is based on the framework proposed 
by Toutain, Elsworth & Chaplin (2006). The influence of 
source localization (Chaplin & Appourchaux 1999) does not 
need to be explicitly considered as it produces a line shape 
that is the same as from the correlated background. We in¬ 
clude it implicitly by our choice of correlation coefficient. 
Thus, there are just two factors in our method that con¬ 
tribute to the asymmetry of the artificial mode peaks. First, 
background noise is correlated with the excitation of the 
modes, and second, overtones of the same angular degree 
and azimuthal order have excitation functions that are cor¬ 
related in time. This correlated mode excitation is based 
on the description given in Chaplin, Elsworth & Toutain 
(2008) and the method is described in detail in the following 
sections. Further details, and analytical descriptions of the 
examples in Sections 2.1 and 2.1 below, may be found in 
Chaplin et al. (2008). 

2.1 Correlated noise and mode excitation 

In Toutain et al. (2006) it was hypothesized that the excita¬ 
tion function of a mode of angular degree I, azimuthal degree 
m, and frequency v, is the same as that component of the so¬ 
lar background (granulation) noise that has the same spher¬ 
ical harmonic projection, in the corresponding range in 
frequency in the Fourier domain. Let us start by considering 
the impact on the line shape of this correlated background 
noise. We do this in the context of how the solarFLAG sim¬ 
ulator generates artificial data on single p-modes. 

The basis of the solarFLAG simulator is the method 
discussed by Chaplin et al. (1997) for generating time series 
of individual p modes. The method uses the Laplace trans¬ 
form solution of the equation of a forced, damped harmonic 
oscillator to make the output velocity of each artificial mode. 
Oscillators are re-excited at each time sample - the chosen 
cadence is typically 40 sec or 60 sec for Sun-as-a-star data - 
with small ‘kicks’ from a time series of random noise. This 
procedure mimics the stochastic excitation of the solar p- 
modes. For the moment we shall assume the random noise 
data are drawn from a normal (i.e., white) distribution hav¬ 
ing zero mean. We shall see later that in the final version 
of the simulator we used ‘granulation-like’ noise - to mimic 
the spectrum of solar granulation - which has a spectrum 
that looks white only locally in the vicinity of the oscillator 
resonance. 

The left-hand panel in Fig. 1 shows, on a logarith¬ 
mic scale, the limit power spectral density of two scenarios, 
which illustrate the impact of correlated noise. In both sce¬ 
narios the power spectral density has contributions from a 
mode - of frequency v = 2990 pHz, and line width F = 1 pHz 
- and from white background noise. The thin dashed red 
line shows the result for a scenario where the excitation of 
the mode is uncorrelated with the background noise. The so¬ 
larFLAG simulator would generate the data for this scenario 
by using one time series of random noise to excite the mode; 
to the velocity output of the mode would then be added an¬ 
other, completely uncorrelated time series of appropriately 
scaled random noise, to represent the background (here the 
height-to-background ratio is 100). Since there is no correla¬ 
tion, the limit frequency power spectrum is given simply by 
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the incoherent addition of the limit spectrum of the mode 
(a Lorentzian) and the limit spectrum of the noise (here, a 
flat offset for white noise). 

Now, what if we use the excitation time series as the 
background noise? The excitation and the background will 
now be 100-per-cent correlated, and the power spectral den¬ 
sity will have a peak that is asymmetric. This situation, for a 
large (about ten per cent) value of the asymmetry, is shown 
with the thick black line in the left-hand panel of Fig. 1. 
If only the effect of the background correlation had been 
considered the peak asymmetry would have been positive 
with an excess of power on the high-frequency side of the 
peak. However, the asymmetry seen in low-Z solar p modes 
observed in Doppler velocity is negative and hence we show 
the effect of negative correlation with excess power on the 
low-frequency side of the resonance. The observed asymme¬ 
try is much lower than illustrated here. 

Next, we consider the impact of correlated mode exci¬ 
tation. 

An important implication of the framework proposed by 
Toutain et al. (2006) is that overtones with the same {I, m) 
should have excitation functions that are correlated in time. 
(Note that the Yi™, for (Z, m) and (Z, —m) are orthogonal, 
and are therefore assumed to have independent, i.e., uncor¬ 
related, excitation.) To illustrate the impact of this correla¬ 
tion, we look at the simplest possible scenarios, which have 
just two modes in the frequency power spectrum. Consec¬ 
utive overtones of the low-Z solar p modes are separated in 
frequency by ~ 135 /rHz. Here, we consider two modes sepa¬ 
rated in frequency by 20 ^Hz. This smaller frequency spacing 
exaggerates the impact of the correlated excitation on the 
observed power spectral density, and therefore allows us to 
show more clearly the effect of the correlation. We again as¬ 
sume each mode has a line width of 1 /rHz. The impact of 
the actual ~ 135-/rHz spacing is considered in Subsection 2.2 
below. 

The thin dashed red line in the middle panel of Fig. 1 
shows the limit spectrum for two modes whose excitation 
is uncorrelated in time. There is no background noise. The 
solarFLAG simulator would generate the data for this sce¬ 
nario by using independent time series of random noise to 
excite each oscillator. Now, what happens if the simulator 
excites both oscillators with the same time series of ran¬ 
dom noise? The excitation is now 100-per-cent correlated 
in time, and the power spectral density shows clearly that 
the peaks are asymmetric (thick dark line). This asymmetry 
comes from the interaction of the tails of the mode peaks. 
We draw an important conclusion from this example: cor¬ 
related excitation of modes will give a contribution to the 
observed asymmetry that is dependent on how the tails of 
the individual modes overlap. It turns out that this effect is 
important even for quite well separated modes. 

In our final example, we consider both sources of asym¬ 
metry together and we add background noise to each two¬ 
mode scenario above. The thin dashed red line in the right- 
hand panel of Fig. 1 is for a scenario where the background 
noise and the mode excitation are all uncorrelated. The thick 
dark line instead shows what happens if we add correlated 
background noise to the two correlated modes. This means 
the excitation and background noise are all 100 per cent cor¬ 
related, and we see that addition of the background further 
modifies the shape of the power spectral density, relative 


to the correlated-mode example with no background noise. 
As such, there are now two factors which contribute to the 
peak asymmetry: there is a contribution from the correlated 
excitation of individual modes (and also potentially from 
source localization); and a contribution from the correlated 
background. 

One other important point to note is that correlation 
of the excitation in time does not imply correlation of the 
mode amplitudes in time. This can be understood by con¬ 
sidering the analogy of damped, stochastically forced oscil¬ 
lators. Modes of different frequencies will be ‘kicked’ by the 
common excitation at different phases in their oscillation cy¬ 
cles, and provided the frequencies differ by more than a few 
line widths (see Chaplin et al. 2008) - a condition easily met 
by consecutive overtones of the low-Z modes, which are sep¬ 
arated by ~ 135 ^Hz - there will be signihcant differences 
in how the amplitudes vary in time, due to the excitation. 

2.2 The solarFLAG simulator 

2.2.1 General overview 

The solarFLAG data sets simulate full-disc ‘Sun-as-a-star’ 
Doppler velocity observations of the Sun, such as those made 
by the BiSON and GOLF. The solarFLAG data sets are 
made with a full cohort of simulated low-Z modes, covering 
the ranges 0 ^ Z ^ 5. The frequencies of the modes come 
from a standard solar model of the user’s choice. A surface 
term is also added to these frequencies, based on polyno¬ 
mial fits to the differences between the standard model fre¬ 
quencies and frequencies from analysis of BiSON and GOLF 
data. For data presented in this paper the frequencies came 
from model BS05(OP) of Bahcall et al. (2005). 

A database of p-mode power and line width and asym¬ 
metry estimates, obtained from analyses of GOLF and Bi¬ 
SON data, was used to guide the choice of the other mode 
input parameters. The hypothetical solarFLAG instrument 
was assumed to make its observations from a location in, 
or close to, the ecliptic plane. This is the perspective from 
which BiSON (ground-based network) and GOLF view the 
Sun. The rotation axis of the star is then always nearly 
perpendicular to the line-of-sight direction, and only a sub¬ 
set of the 2Z + 1 components of the non-radial modes are 
clearly visible: those having even Z -|- m. These components 
are represented explicitly in the solarFLAG time series. 
The visibility for given (Z, m) also depends, although to a 
lesser extent, on the spatial filter of the instrument (e.g. 
Ghristensen-Dalsgaard 1989; Broomhall et al. 2009). Here, 
we adopted BiSON-like visibility ratios. 

Fig. 2 shows a schematic representation of the so¬ 
larFLAG simulator, the details of which we discuss next. 

Before the simulator can begin its calculations, all the 
input parameters must be specified. The p-mode input pa¬ 
rameters are held in a control file (with mode powers having 
been suitably scaled to reflect the visibility filter of the Sun- 
as-a-star observations). The other input parameters relate 
to the granulation-like noise, and other background noise, 
and are specified at run time by command-line inputs. Time 
series of granulation-like kicks are used to excite the modes, 
and are also used to give correlated background noise. The 
timescale and standard deviation of this noise must be cho¬ 
sen. A single constant, p, is also set and this fixes the co- 
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Figure 1. Limit frequency power spectra for scenarios with one or two modes. Left-hand panel: a single mode, with background noise. 
The dashed red line shows the spectrum for no correlation; the solid black line tor when the mode is correlated with the background 
noise. Centre panel: two modes. The dashed red line shows the limit spectrum when the excitation of the modes is uncorrelated; the 
solid black line when the excitation is correlated. Right-hand panel: Two modes, with background noise. The dashed red line shows the 
limit spectrum when there is no correlation of the excitation, or with the background noise; the solid black line when the excitation is 
correlated, and there is correlation with the background noise. 


efficient of correlation for the correlated excitation, and the 
correlation with the background noise. This gives the user 
the flexibility to ‘tune’ the asymmetry of the mode peaks 
(as well as other effects arising from the correlations), since 
higher correlation gives larger peak asymmetry. The input 
parameters of other sources of noise, such as the standard 
deviation, CTpsn, of the photon-shot noise, are also specihed 
at this point. 

The main block of the solarFLAG code generates out¬ 
put time series for each (Z,m). It makes the time series of 
all the overtones, n, for that {l,m). The modes are excited 
by time series of granulation-like noise. For each [I, m) there 
is a noise time series that gets used again and again to ex¬ 
cite the overtones: this is the correlated gr. noise time series. 
When the user fixes the correlation at less than 100 per cent, 
a suitable fraction of uneorrelated gr. noise (see schematic 
in Fig. 2) must be mixed in when the modes are excited. 
These uncorrelated noise time series must be seeded with 
completely independent random noise and made afresh for 
each mode. After all overtones have been made, a suitably 
scaled mixture of correlated and uncorrelated noise is added 
to the time series. This gives the cumulative time series for 
each Execution of the main code block is then re¬ 

peated for a total of 14 different {I, m) combinations^. 

After time series for all the {I, m) have been made they 
are added together. The final stage of processing adds in the 
simulated time series of other uncorrelated noise sources. 
The schematic includes a contribution from photon-shot 
noise. Other sources may also be included at this point (e.g., 
active-region noise and instrumental noise) but this was not 
done for the data used in this work. 

A complete description of the different elements that 
make up the complex frequency amplitude (and frequency 
power) spectrum of a solarFLAG time series is presented in 
Appendix A. 


^ Data on modes cover the range 0 Si Z Si 5; and there are I 1 
m-components visible in the Sun-as-a-star data at each 1. This 
actually sums to give 21 combinations of However, com¬ 

ponents at Z = 4 and 5 that do not have Z = ±m are so weak that 
they are not included in the time series generation. Discounting 
these modes reduces the number of combinations to 14. 


We now turn to specifics on how the data are made, and 
how the effects of correlations are quantified. 


2.2.2 Achieving correlated excitation of overtones 

As noted above, we use time series of granulation-like noise 
to excite the overtones, n, of the same I and m. Later we 
will refer to this as ‘correlated noise’. The granulation-like 
noise is made using white-noise input to a low-order, au¬ 
toregressive model of the AR[1] type (e.g., see Koen 2005), 
i.e.. 


u{t) = u{t — At)exp—At/T + S{t). (1) 


Here, u{t) is the output time series of granulation-like noise; 
T is the time constant of the model, which should be given 
a value to mimic the lifetime of the solar granulation (see 
below); At is the cadence at which samples are generated 
by the autoregressive process; and finally the S{t) are ran¬ 
dom numbers drawn from a normal distri bution h aving zero 
mean and sample standard deviation ybr^At/r, where a 
fixes the amplitude of the granulation-like noise. The power 
spectral density, n(iz), of this granulation-like noise follows 
the approximate but useful Harvey (1985) power-law model 
i.e.. 


n{iy) 


4a^T 

1 -|- {2'kvtV 


( 2 ) 


The default option would be to excite all overtones of the 
same {l,m) with the same time series u{t), meaning the ex¬ 
citation of the overtones would be 100-per-cent correlated 
(as in Section 2.1). As indicated earlier, we decided to add 
to the solarFLAG simulator the option to let the user choose 
a coefficient of correlation, p, common to all modes. We de¬ 
scribe in Section 2.2.2 below how the asymmetry of a mode 
is fixed by the combination of p, the noise background, and 
the input mode paraeters. 

Overtones of the same (I, m) are excited by a composite 
time series of kicks, u{t), where: 

w(i) = |p|^^^“c(t) + (1 - \p\f^^ Unit). (3) 


In the above, the Uc{t) and Uu{t) are both time series of 
granulation-like noise. However, the Uc{t) are kicks that are 
common to all the overtones (the correlated gr. noise time 
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Figure 2. A schematic representation of the solarFLAG simulator. 
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series shown in red in Fig. 2); while the Wu(t) are completely 
independent (shown as uncorrelated gr. noise in Fig. 2). 

So far in this section we have not said anything about 
tuning the correlation with the background noise. Correlated 
background noise is provided by the correlated part of the 
u{t). After all overtones of the same {l,m) have been made, 
the u{t) is added to the time series to give background noise 
(having first been appropriately scaled ). In the formulation 
presented by Chaplin et al. (2008), another coefficient, a, 
was used to describe the correlation with the background 
noise. Here, we assume this correlation has the same size as 
the correlation between the excitation of different modes, so 
that we set a = p. The a priori choice of p therefore also 
hxes the correlation of the overtones with the background 
noise. 

How do we juggle all the choices implied by the descrip¬ 
tions above to give a time series where both the S/N in the 
modes and the asymmetry of the modes are realistic? That 
is the question we turn to next. We use two main factors to 
contribute to the peak asymmetry: correlation of the modes 
with the background noise; and correlated excitation of the 
overtones. We will now consider each of these in turn as it 
relates to controlling the asymmetry of the artificially gen¬ 
erated modes. 


2.2.3 Asymmetry due to correlated noise 

The contribution to the asymmetry from correlations with 
the noise background may be dealt with analytically in a 
fairly straightforward manner. Let fenim be the asymmetry, 
due to correlated noise, of the nth overtone of a given (Z, m). 
The frequency of this mode is Onim. The asymmetry (as 
it appears in the mode parametrization function in Equa¬ 
tion 12) is given by (Toutain et al. 2006): 


nimiynl D 


H„i n 


(4) 


where Hnim is the height of the mode peak in the frequency 
power spectrum. Notice the presence of the factor p in Equa¬ 
tion 4: it is needed because it is only the correlated part of 
nimionim) that contributes to the asymmetry. nimVnim) is 
the granulation-noise background at the frequency of the 
resonance, which is just (cf. Equation 2): 


r^lmVnlm) — 


1 -I- {2%Vnl rn'T 


(5) 


We give aim in Equation 5 an explicit dependence on the 
combination {l,m). This gives us the ability to choose dif¬ 
ferent aim for different (Z, m), to take proper account of the 
different mode visibilities. This plays a key part in tuning 
the asymmetry, as we now go on to explain. 

We have assumed that at a given frequency the rela¬ 
tive sizes of the granulation noise and the mode amplitudes 
on the Sun are independent of degrees Z and m. An impor¬ 
tant consequence of this assumption is that the asymmetries 
from the correlated noise will be the same for all (Z,m). In 
the solarFLAG simulator, after making overtones of a given 
(Z,m), the time series of mode amplitudes is multiplied by 
a visibility factor, Sim, to mimic the visibility filter of the 
Sun-as-a-star observations. To preserve the independence of 
asymmetry on (Z, m) we must therefore also modify the aim 


of the granulation noise, according to: 

Sim 


^lr\ 


Sim 

all {l,m) 


(6) 


Here, a is the equivalent standard deviation of the cumula¬ 
tive granulation noise background, n{v), so that: 


i{u) = 


- E 


^afmT 


1 -I- (2'kvtV 1 + (2'kvtV 

all {l,m) 


(7) 


We have now covered all the steps needed to see how the 
asymmetry is tuned. 

The p-mode parameters (i.e., frequencies, heights, line 
widths) are fully specified on input. The asymmetry is then 
tuned by three free parameters: the coefficient of correlation 
p; the equivalent standard deviation of the cumulative gran¬ 
ulation noise background, a\ and the timescale of the granu¬ 
lation, r. Once a and r have been chosen, we use Equation 6 
to determine the aim for each combination of (Z, m). This in 
turn determines the asymmetry, due to the correlated noise, 
according to Equations 4 and 5. 

The left-hand panel of Fig. 3 shows resulting asymme¬ 
tries given by the parameters a = 0.2 m s“^ and r = 260 sec. 
The dotted line shows results for p = —0.20, the solid line 
for p = —0.36, and the dashed line for p = —1.00. Here, we 
used negative p to give the negative asymmetries observed 
in Doppler velocity data. The magnitude of the asymmetry 
increases in direct proportion to the magnitude of p is in¬ 
creased (this follows trivially from inspection of Equation 4) . 
When p = —0.36 we get asymmetries that resemble quite 
closely those seen in the real observations. The asymmetries 
are largest at low and high frequency, because the ratio of 
the granulation background to mode amplitude is highest 
there. Very similar looking plots are given by varying cr; in¬ 
spection of Equations 4 and 5 indicates that the magnitude 
of the asymmetry is proportional to a. 

The left-hand panel of Fig. 4 shows instead the effect 
of varying r (for fixed a = 0.2ms“^ and p = —0.36). The 
dotted line shows results for r = 130 sec, the solid line for 
T = 260 sec and the dashed line for r = 520 sec. Changes in 
T clearly have a less pronounced effect on the asymmetries 
than do similar relative changes in p. 


2.2.Asymmetry due to correlated excitation 

The contribution to the asymmetry overtones of the same 
{l,m) from the overlapping wings of the overtones of the 
same (Z, m) that have been excited with correlated noise is 
fixed by the frequency separations, line widths and relative 
heights of those overtones, and the choice of p. We later 
refer to this as ‘correlated excitation.’ This contribution is 
far less amenable to a neat analytical description than the 
correlated noise contribution described above. For a given 
overtone it must describe, in the complex plane, the con¬ 
tribution of power from all other overtones to which it is 
correlated. A full description of the complex frequency am¬ 
plitude spectrum, and frequency power spectrum, is given in 
Appendix A. We may use the analytical descriptions there 
to calculate numerically the asymmetries given to the mode 
peaks by the correlated excitation. The resulting estimates 
are shown in the central panels of Fig. 3 and 4. Note that 
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From: correlated noise 


From: correlated excitation 


Total 



Figure 3. Impact of choice of p on mode peak asymmetry. The left-hand panel shows the asymmetry due to correlation with background 
noise (Section 2.2.3); the centre panel the asymmetry due to correlated mode excitation (Section 2.2.4); and the right-hand panel the 
total asymmetry, due to the complex interaction of the two contributions. All panels show results for a = 0.2 ms“^ and r = 260 sec. 
Different linestyles show results for different p: dotted for p = —0.20; solid for p = —0.36 and dashed for p = —1.00. 


From: correlated noise From: correlated excitation Total 



Figure 4. Impact of choice of r on mode peak asymmetry. The left-hand panel shows the asymmetry due to correlation with background 
noise (Section 2.2.3); the centre panel the asymmetry due to correlated mode excitation (Section 2.2.4); and the right-hand panel the 
total asymmetry due to the complex interaction of the two contributions. All panels show results for a = 0.2ms“^ and p = —0.36. 
Different linestyles show results for different r: dotted for r = 130 sec; solid for t = 260 sec and dashed for r = 520 sec. 


the correlated excitation contribution does not depend on 
the sign of p (only the magnitude), and is independent of 
the choice of r (and a). 

At lower frequencies the asymmetry given to modes by 
the correlated excitation is negative. Here, there are a larger 
number of other overtones at frequencies above a mode than 
there are below it. The reverse is true at higher frequencies, 
where the asymmetry is positive. To understand the sign of 
the asymmetry in each region, we refer back to the central 
panel of Fig. 1. There, we had two correlated modes. The 
lower-frequency mode had negative asymmetry, due to the 
correlated impact of its higher-frequency counterpart; the 
higher-frequency mode showed asymmetry of opposite sign, 
due to the lower-frequency mode. We see this pattern re¬ 
peated in the full spectrum of overtones with the cross-over 
in behaviour located where p-mode power is a maximum. 

The right-hand panels of Fig. 3 and 4 show the full 
asymmetry given to the modes. It is worth pointing out 
that this is not simply the sum of the asymmetries given 
by correlated noise (left-hand panels) and correlated excita¬ 
tion (central panels). Matters are a bit more complicated, 
because the correlated noise and correlated excitation act 
together in a non-trivial manner in the complex plane. 

One final thing to mention with regard to the asymme¬ 
tries is that there is actually a third contribution, which 
comes from the fact that the frequency response of the 
granulation-like noise used to excite modes is not white 
(i.e., flat). The response in the vicinity of each resonance of 
course rises with decreasing frequency, meaning there will 


be a small negative asymmetry contribution. The effect is, 
however, negligible. (The impact of the non-white response 
of the excitation is modelled in Appendix A.) 

We have seen that the input frequencies, powers and 
line widths of overtones, and the correlation coefficient p, 
fix the asymmetry contribution from the correlated excita¬ 
tion. When in addition the cr and r of the granulation-like 
noise are chosen the asymmetry contribution from the corre¬ 
lated noise is also completely specified. Once any additional 
uncorrelated background has been specified, these choices 
in principle also fix the observed S/N ratio of the mode 
peaks. For the data analysed below, we chose r = 260 sec 
- the characteristic solar timescale - and p — 0.36, chosen 
to give a good match to the asymmetry measured in pre¬ 
vious work. These values correspond to the solid curves in 
Figures 3 and 4. We also chose a = 0.2 ms“^ for corre¬ 
lated and a = 0.25 ms“^ for uncorrelated noise. We now 
go on to describe the S/N, in terms of the commonly used 
background-to-height ratio, /3. 


2.2.5 Background-to-height ratio in mode peaks 


The intrinsic background-to-height ratio, pnim, in a mode 
peak is given by: 


CnlTn) 

Hr,l m 


(8) 


where N{vriim) is the total background at the resonance. It 
has a part due to the cumulative granulation background of 
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all [I, m) that appear in the time series. There will also be 
parts due to other sources of uncorrelated noise. In what fol¬ 
lows we specify just one of the possible components: photon 
shot noise, A^pan. This shot noise is specified by its variance, 
f’'psn- The power spectral density due to this noise term is 
then: 

iVpan =2a = anAt. (9) 

The total background may therefore be written: 

N {ynlm) ~ A^psn T ^ ^ -t (^0) 

all {l^rri) 

or explicitly: 

= 2a^an At + ^ ^ ■ (11) 

The above implies that the intrinsic background-to-height 
ratio is in principle specified by the choice of four param¬ 
eters (when no instrumental or other sources of noise are 
specihed): the height of the mode, Hnim', the a and r of 
the granulation-like noise background; and the (Jpsn of the 
uncorrelated shot-noise background. 

2.2.6 Solar-cycle effects 

It is a fairly straightforward matter to include solar-cycle¬ 
like variations in the artihcial timeseries data. Changes in 
frequency may be introduced by varying in time the natu¬ 
ral frequency of the damped, driven oscillator used to sim¬ 
ulate each mode component; whilst changes in amplitude 
and damping may be introduced by varying in time the 
damping constant for the oscillator. In all cases we require 
changes that mimic those seen in the real data. We base 
these changes on an artificial proxy of the real 10.7-cm radio 
flux variations observed over solar activity cycle 23. Varia¬ 
tions in the frequencies and damping rates are programmed 
to follow the artihcial proxy, with appropriate calibration 
constants used to hx the absolute scale of the variation for 
each parameter. 

Here, we used a high-order polynomial ht to the 10.7- 
cm radio flux versus time - for the 11-year period beginning 
1997 February 3 - to provide the proxy activity, which we 
refer to as X in later sections of the paper. The purpose 
of the fit was to provide a smooth proxy that captures the 
main 11-year cycle while removing shorter-term variations. 
Calibration constants for the frequencies and damping rates 
of every simulated mode were based on results from previous 
analyses of real BiSON data. The calibration constants for 
the frequencies depend not only on frequency - the higher 
the frequency, the larger is the solar-cycle frequency shift 
- but also the angular degree, I, and the azimuthal order, 
m, of each simulated component (e.g. Chaplin et al. 2004). 
The dependence on {l,m) is due to the strong latitudinal 
dependence of the real near-surface activity, which drives 
the changes in the mode parameters. 

As noted earlier, due to the relatively large uncertain¬ 
ties in the measurements little information is available on 
the frequency and {l,m) dependence of detected solar-cycle 
changes in the mode amplitudes and damping rates for low- 
degree modes, but we know from medium- and high-degree 
studies (e.g. Komm et al. 2000a; Howe et al. 2004) that the 


changes are strongest in the middle of the five-minute band. 
However, the relative sizes of the variations are consistent 
with the hypothesis that both are the result of net changes to 
the damping only (e.g. Chaplin et al. 2000). Here, we have 
introduced solar-cycle variations that adhere to this finding, 
so that in our simulations we change only the damping con¬ 
stant. Further information on the introduced variations is 
included later in Section 4. 

Finally, we note that variations in amplitude lead to 
solar-cycle-like changes in the mode peak asymmetries of 
our artificial data. This is because the height-to-background 
ratios of the mode peaks changes with the artificial proxy, 
due to variations in the intrinsic amplitudes (the simulated 
granulation and shot-noise components were stationary in 
time). Again, information on the programmed variations is 
in Section 4. 


3 FITTING METHODS 

Part of the object of the current work is to compare the 
results from two different fitting algorithms. Both of these 
are conventional in the sense that they use the “pairwise” 
fitting approach in which each 1 = 0/2 and 1 = 1/3 pair of 
peaks is fitted independently with its own background and 
parameter set. 

The two methods we use are labelled as MLE (Max¬ 
imum Likelihood Estimate) and MCMC (Markov Chain 
Monte Carlo) but are not solely described by the method of 
optimization. There are differences in the power spectrum 
model that we detail below. 

3.1 Mode Parametrization 

Common to both fitting methods, each mode component of 
degree I, radial order n and azimuthal order m is represented 
by a peak profile described by the equation 

^(^) = (r^) ^ [(1+''^)'!’ (12) 

where, 

c = 2 ( 1 / - oo)/r, (13) 

i/Q is the frequency of the Lorentzian component, F its 
width, h its height, and b a fractional parameter char¬ 
acterizing the asymmetry. This expression simplifies to 
the normal Lorentzian for b = 0. This is the profile of 
Nigam & Kosovichev (1998), with the quadratic term in b 
suppressed for greater stability (Fletcher et al. 2009). For 
this work, we have chosen to parametrize the peak in terms 
of the amplitude. A, where 

h = A^/lnV). (14) 

Neglecting the small asymmetry term, JV = ttF/i is propor¬ 
tional to the integrated energy of the mode (see, for example, 
Komm et al. 2000b, but note that their A is our h), and the 
energy supply rate dE/dt is given by 

dE/dt oc A^T oc Vh. (15) 

The likelihood spaces, (or the posterior probability dis¬ 
tributions) are log-normal for the mode width and mode 
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Figure 5. Sensitivity of fitted frequency to activity level and fill for artificial data. Panels (a)—(d) show fitted (symbols) and input (lines) 
dUnildX for MLE (a, c) and MCMC (b,d) fits to 11 one-year sets of artificial data at different activity levels, mimicking the solar cycle. 
The results are for 100 per cent (a, b) and realistic (c, d) duty cycle. The lower two rows, for MLE only, show the results of the tests 
with multiple realisations of the artificial data. Panel (e) shows the sensitivity of the frequency to activity index, X\ panel (f) is similar 
to panel (e) but for a set of artificial time series where there is no input frequency variation. Panel (g) shows the sensitivity to duty cycle, 
E, for the test with 25 realisations of the artificial data. Panel (h) shows the rate of change with duty cycle F of the sensitivity dvjdX 
of mode frequency to activity index, for the test where each of the 21 one-year BiSON window functions was applied to each of the 11 
years of the SolarFLAG time series. Black circles represent / = 0, blue diamonds / = 1, red squares 1 = 2, and green open triangles 1 = 3. 
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amplitude parameters, so we follow the common practice of 
varying the logarithms of the amplitude and width parame¬ 
ters rather than their raw values. 

The background term is a constant for each Z = 0/2 or 
I = 1/Z pair, and in both algorithms the asymmetry term a 
is also common to all of the peaks in each pair. In the code 
used for the MLE algorithm there is a separate T for each n 
and I, while the MCMC code uses the same T for every peak 
within anZ = 0/2orZ = l/3 group. There is a separate am¬ 
plitude parameter for each n and I, but the relative heights 
of the rotationally split components of different m within 
an n,l group are fixed in the MLE code at 1:0.41 for the 
m = 2 : m — 0 component ratio for I = 2 and 1 : 0.19 for the 
m = 3 : m = 1 ratio for 1 = 3. In the MCMC code they are 
allowed to vary; the priors are set at a normal distribution of 
width 0.2 centered on 0.4 for the I — 2,m = 0 ■. I = 2,m = 2 
ratio and the I — 3, m = l-. I — 3, m = 3 ratio while for 
the I — 3, m = l: I — 3, m = 3 ratio we use a flat-topped 
Gaussian centered on 0.2 that is flat between 0.1 and 0.3 
and has Gaussian sides that drop off with width 0.1. As we 
will see, the differences have no obvious impact on the final 
results. 

The rotational splitting parameter was fixed at 0.4 pHa 
for the MLE algorithm, while for the MCMC it was left as 
a free parameter with a Gaussian prior with centre 0.4 and 
width O.OS/rHz. 


3.2 Optimization 

The Fourier power spectrum has the statistics of with 
two degrees of freedom rather than Gaussian statistics. The 
quantity to be mimimized is the negative log of the likeli¬ 
hood function (Anderson et al. 1990), 

5 = -lnL = ^{lnMi + ^}, (16) 

i 

where M is the model, O is the observed spectrum, and the 
sum is over the data points (frequency bins) in the fitting 
window. 

In the MLE approach we simply find the values of the 
model parameters that give the minimum value of S and 
estimate the errors by inverting the Hessian matrix. For 
MCMC, on the other hand, we take the standard deviation 
of the posterior distribution for each parameter. 

The Maximum Likelihood Estimation (MLE) code used 
here was a slightly modified version of the one developed 
for use with BiSON data and described by Chaplin et al. 
(1999). The MCMC code is a vanilla Markov Chain Monte 
Carlo approach (see Gelman et al. 2003; Davies et al. 2014a, 
for more details). 


3.3 Treatment of window function 

It is a well-known problem in ground-based helioseismology 
that the Earth’s rotation causes the observations to be inter¬ 
rupted with a 24-hour periodicity, which gives rise to ‘side- 
lobe’ peaks at multiples of 11.57 pHz from each mode. At 
low degrees, where I = Q/2 and Z = 1/3 pairs are separated 
in frequency by a similar amount, this is particularly prob¬ 
lematic. The effect is greatly reduced by the use of multi-site 


observations but cannot be entirely eliminated, and it must 
therefore be taken into account in the data analysis. 

As is conventional, we describe the modulation of the 
time series by a ‘window function’ - a time series corre¬ 
sponding to the observations in which each time sample is 
represented by a one for an observation and a zero for miss¬ 
ing data. The observed spectrum can be considered as the 
convolution of the power spectrum of this window function 
with that of the uninterrupted observations. 

In the MLE code used here, the sidelobe peaks are rep¬ 
resented in the htting model by peaks of the same width and 
asymmetry as each main peak, at ±11.57 pHz from its cen¬ 
tral frequency. The initial guess for the relative amplitude of 
this sidelobe peak was estimated from a spectrum in which 
a single sine wave was convolved with the window function 
of the time series. For modes in the main five-minute band 
(between 2.0 and 3.5 mHz) this relative height was then al¬ 
lowed to vary in the first pass of htting and held Hxed for 
modes outside this frequency range. The mean value of the 
htted sidelobe after the hrst pass was then used as the hxed 
value for all modes in the hnal pass. This allowed a larger 
number of successful hts to be recovered than if the sidelobe 
height had been left as a free parameter throughout. 

The MGMC algorithm uses a different, more accurate 
but more computationally expensive approach in which the 
model spectrum is convolved with the window function spec¬ 
trum at each optimization step. This method can also ac¬ 
count for the broadening of the peak as the duty cycle of 
the observations is decreased; however, neither this approach 
nor the one used with the MLE algorithm take into account 
the correlations between frequency bins introduced by the 
missing data. Note that the treatment of the sidelobes is 
related to the specifics of the codes used here and is not 
intrinsic to the MLE or MCMG approach; in principle it 
would be possible to perform MLE with the convolution, as 
was indeed done for example by Jimenez-Reyes et al. (2007); 
Fletcher et al. (2009). 


4 ARTIFICIAL DATA TESTS 

The object of the tests described here is to check that the 
MLE and MCMC algorithms give consistent and reliable 
results for the variation of the different mode parameters 
with the activity index. 

The 11-year time series generated as described in Sec¬ 
tion 2 above was divided into contiguous 365-day segments 
and each segment was fitted using both the MLE and 
MCMC algorithms. The series were fitted both with 100 
per centduty cycle and with a ‘realistic’ duty cycle based on 
that of the BiSON time series for the 11 years starting with 
the notional start date of the artificial series, 1997 Febru¬ 
ary 3. For each mode we then take an error-weighted mean 
over all the data sets and subtract this from the values to 
obtain the shifts, finally using a linear least-squares ht to 
the activity index to derive the sensitivity for each mode for 
comparison with the input values. 

For the MLE method only, three additional, more ex¬ 
tensive tests were carried out. In one test, 25 additional re¬ 
alizations of the 11-year artificial data series were fitted in 
one-year segments, while in another each one-year segment 
of a single realization of the artificial data was fitted with the 


MNRAS 000, 000-000 (0000) 


12 R. Howe, G.R. Davies, W.J. Chaplin, Y.P. Elsworth, and S.J. Hale 




c d 







h 



1500 2500 3500 4500 

L' (/zHz) 


Figure 6. Sensitivity of mode amplitude to activity index and duty cycle. Panels (a)-{d) show fitted (symbols) and input (lines) 
d\ogAni/dX for MLE (a, c) and MCMC (b,d) fits to 11 one-year sets of artificial data at different activity levels, mimicking the solar 
cycle. The results are for 100 per cent (a, b) and realistic (c, d) duty cycle. The bottom two rows, for MLE only, show the results of the 
tests with multiple realizations of the artificial data. Panel (e) and (f) show the sensitivity of the amplitude to activity index, X, with 
panel (f) showing the results for the test with no input frequency variation. Panel (g) shows the sensitivity to duty cycle, F, for the test 
with 25 realizations of the artificial data. Panel (h) shows the rate of change with duty cycle F of the sensitivity d log A^i/dX of mode 
amplitude to activity index, for the test where each of 21 one-year BiSON window functions was applied to each of the 11 years of the 
SolarFLAG time series. Black circles represent I = 0, blue diamonds / = 1, red squares 1 = 2, and green open triangles I = 3. 
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Figure 7. Fitted amplitude (left) and linewidth (right) shifts averaged over 25 realizations of the artificial data with realistic duty cycle 
for each mode between 2.0 and 3.5 mHz, plotted against the input value. Black circles represent 1 = 0, blue diamonds 1 = 1, and red 
squares 1 = 2. The correspondingly coloured lines represent least-squares fits for each 1. 


duty cycle corresponding to each of the 21 one-year BiSON 
time series. The term ‘realization’ here refers to a realiza¬ 
tion of the noise exciting the simulated modes. The extra hts 
make it easier to distinguish the systematic errors from the 
random noise; by using 25 realizations we reduce the ran¬ 
dom errors by a factor of five, which is sufficient to clarify 
the trends. In the real observations we have 22 years of data, 
so the random errors on the mean parameters and sensitiv¬ 
ities are a factor of \/2 smaller than for a single realization 
of the simulated data, or approximately three times larger 
than for 25 realizations. Finally, to help rule out the pos¬ 
sibility that the frequency shifts are caused by asymmetry 
changes, we prepared another set of artificial data realiza¬ 
tions where the input amplitude and width changes - and 
hence the asymmetry changes - were as before but the input 
frequency change was set to zero. We note that to generate 
frequency shifts as large as those observed, the asymmetry 
shift would need to be extremely and unrealistically large, 
as discussed below in Section 4.2. 


4.1 Results 

4.1.1 Frequencies 

Figure 5 (a-d) shows the input and fitted values of dvni/dX, 
the frequency shift per unit of the activity index X for 
each fitting method. Note that the differential notation we 
use throughout is a convenient shorthand for the regression 
slope. In the case of the artificial data the linear relation¬ 
ship between activity and mode parameters is designed to 
be exact; the real data may exhibit small deviations from 
this but it is still a useful measure of the sensitivity to 
first order. Except for a few modes at the higher end of 
the frequency range, the results from both of the algorithms 
agree within errors with the input values; the mean ratio 
between the fitted and input dv^i/dX for 0 ^ I ^ 2 and 


2000/rHz ^ i/ni ^ 3500/rHz is consistent with unity in each 
of the four cases. 

Panels (e) and (f) of Figure 5 illustrate the sensitivity 
of the frequency to activity for the tests with 25 realizations 
of the artificial data, with panel (f) showing the results for 
the test with no input frequency variation. Apart from a 
few modes at frequencies above 3.7 mHz where the s/n ra¬ 
tio is low, we see no frequency change in the dataset where 
none was explicitly introduced. Figure 5(g) shows the sen¬ 
sitivity of the frequency to the duty cycle. These results 
were obtained by performing for each mode a multiple lin¬ 
ear regression in which the independent variables were the 
activity index X and the duty cycle F and the dependent 
variable was the mode frequency. The results in Figure 5(e) 
clearly show the slightly lower dv/dX for the I = 0 modes - 
a difference that is lost in the noise for the single-realization 
test. 

In another test of the MLE fitting we applied each of 
the 21 one-year window functions from the BiSON network 
for 1993-2013 to each of the 11 years of the SolarFLAG 
time series. These results were used to check for any bias 
in the sensitivity of the frequency to activity level due to 
the different window function, as follows. The results were 
divided into 21 sets in which the 11 one-year spectra all 
had the same duty cycle but different activity indices. For 
each set a regression was performed for each mode, with the 
frequency shift as the dependent and the activity index as 
the independent variable, to obtain a set of slopes dv^i/dX 
for each value of F. Finally, a regression was performed for 
each mode with the duty cycle F as the independent variable 
and the dvjdX value as the dependent one. The results are 
shown in Figure 5(h). The sensitivity of the frequencies to 
activity index shows no systematic bias due to the window 
function. 
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Figure 8. Sensitivity of fitted line width to activity index and fill, for artificial data. Panels (a)-(d) show fitted (symbols) and input 
(lines) dlogFni/dX for MLE (a, c) and MCMC (b,d) fits to 11 one-year sets of artificial data at different activity levels, mimicking the 
solar cycle. The results are for 100 per cent (a, b) and realistic (c, d) duty cycle. The bottom two rows, for MLE only, show the results 
of the tests with multiple realizations of the artificial data. Panels (e) and (f) show the sensitivity of the linewidth to activity index, X, 
with the results in panel (f) being for the test with no input frequency vrition. Panel (g) shows the sensitivity to duty cycle, F, for the 
test with 25 realizations of the artificial data. Panel (h) shows the rate of change with duty cycle F of the sensitivity dF^i/dX of mode 
width to activity index, for the test where each of the 21 one-year BiSON window functions was applied to each of the 11 years of the 
SolarFLAG time series. Black circles represent / = 0, blue diamonds I = 1, red squares 1 = 2, and green open triangles I = 3; for MCMC 
black circles represent / = 0/2 pairs and blue diamonds I = 1/3. 
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4 .. 1.2 Amplitudes 

For mode amplitude A (defined as in Equation 14 above) 
and width F we work with the natural logarithm of the pa¬ 
rameter; a shift in the logarithmic parameter is equivalent 
to a fractional shift in the raw parameter. 

Figure 6(a-d) shows the sensitivity, dlog A^i / dX , of the 
natural logarithm of the amplitude parameter Ani to activ¬ 
ity index, obtained using the two methods for artificial data 
with 100 per cent and realistic duty cycles, compared with 
the input values. The individual dlogAni/dX for the two 
methods are well correlated, with correlation coefficients of 
0.77 for the 100 per cent duty cycle and 0.70 for the realistic 
one over the 47 modes common to both sets; the threshhold 
for 0.1 per cent probability of the same result arising for two 
unrelated populations of 47 samples is 0.47. 

The mean ratios between fitted and input dA„i / dX val¬ 
ues are 0.84 ± 0.08 for MCMC and 0.85 ± 0.07 for MLE in 
the 100 per cent duty cycle case and 0.83 ± 0.09 for MCMC 
and 0.79 ± 0.10 for MLE in the realistic duty cycle case. 
This result suggests that the sensitivity of the amplitude to 
the activity index is being systematically underestimated by 
both methods. However, as we shall see below, when we av¬ 
erage over all realizations of the artificial data we find that 
the average bias is not as severe as for this case. 

In panels (e) and (f) of Figure 6 we show the sensitiv¬ 
ity of the amplitude measurements to activity index for the 
multiple-realization tests with MLE fitting; panel (f) shows 
the results for the artificial data with no input frequency 
variation. The lack of frequency variation does not signifi¬ 
cantly affect the amplitude variation. The results with mul¬ 
tiple realizations clearly reveal the frequency-dependence of 
the shifts, but also confirm the systematic underestimation 
of the sensitivity. The sensitivity to duty cycle (panel (g)) 
has a small non-zero value that does not depend strongly 
on frequency; this is most likely due to the power in higher- 
order daily sidelobes that is not accounted for in the fitting. 

Figure 6(h) shows the variation with duty cycle of the 
amplitude shift per unit activity, from the test with MLE 
fitting where the window function was varied. Apart from a 
few low-frequency modes and the I = 3 modes, the results 
show no systematic bias of the sensitivity measurement due 
to the window function. 

The apparent underestimate of the amplitude shifts is 
of some concern. However, when we consider the MLE fits 
to the larger cohort of artificial time series, we obtain a 
more favourable result. Figure 7(a) shows the fitted ampli¬ 
tude shifts for data with realistic fill, averaged over 25 re¬ 
alizations of the artificial data for each mode and plotted 
as a function of the expected shift. In this case we find, for 
modes between 2.0 and 3.5 mHz, that the slope of a linear 
fit of fitted vs. input shifts is 0.92 ± 0.03, 0.95 ± 0.02, and 
0.89 ± 0.02 for I — 0, I = 1, and I = 2, respectively, which 
suggests that the realization chosen for the main plots was a 
particularly ‘unfortunate’ example and the real systematic 
underestimate of the sensitivity is likely to be less than ten 
per cent. We emphasise that the size of the fractional shifts 
is not affected by any systematic errors in the underlying 
values as long as these are not activity-dependent. 


4.1.3 Line width 

Figure 8(a-d) shows the input and fitted values of 
d log / dX, the rate of change of line width with the activ¬ 
ity index X, for the two methods. As for the amplitude mea¬ 
surements, the individual-mode sensitivity results for the 
two methods are well correlated with one another, with a 
correlation coefficient of 0.59 in the 100 per cent duty cy¬ 
cle case and 0.71 for the realistic duty cycle, over 27 values. 
The slope of a fit of MCMC vs MLE values is 0.96 ± 0.48 
for 100 per cent duty cycle and 1.08 ±0.70 for realistic duty 
cycle, which indicates that the values are consistent; fur¬ 
thermore, for each method the dF/dX values for the 100 
per cent and realistic duty cycle cases are consistent. The 
results of the tests with multiple realizations are shown in 
Fig. 8(e,f). Again, the frequency variation of the sensitivity 
is clearly seen in Fig. 8(e) and (f), with panel (f) being for 
the case with no input frequency variation. Again, the lack 
of frequency variation has no significant effect on the width 
variations. Fig. 8(g) shows that there is a small, more or 
less frequency-independent, effect of the duty cycle on the 
fitted line width, which is to be expected for the MLE al¬ 
gorithm due to the crude handling of the window function 
effects. Fig. 8(h) shows that the results of the MLE test with 
multiple duty cycles applied to the same realization of the 
artificial data; the duty cycle has no significant effect on the 
sensitivity of the linewidth to activity-related variations, so 
as long as the two factors are not correlated we can treat 
the effects as independent of one another. 

As in the case of the amplitude, dF/dX appears sys¬ 
tematically underestimated, with a mean ratio between fit¬ 
ted and input values of 0.81 ± 0.09 for 100 per cent and 
0.79T0.10 for realistic duty cycle for MLE fits and 0.84T0.08 
for 100 per cent fill and 0.83 ± 0.09 for realistic duty cycle 
from MCMC. For the test with MLE fits to 25 realizations 
of artificial data with realistic duty cycle (Fig. 7(b)), we ob¬ 
tain slopes of 0.87 ± 0.04, 0.88 ± 0.03, and 0.77 ± 0.03 for 
I = 0, I = 1 and I — 2 respectively; as with the amplitude, 
this suggests that our main sample realization was particu¬ 
larly unfavourable, but we might still be failing to recover 
the I = 2 shifts completely. 

4 . 1.4 Energy supply rate 

We recall from Equation 15 that the energy of a mode is 
proportional to the square of Its amplitude, and the rate of 
energy supply to the mode dEni/dt is proportional to h„iF^i 
so for the logarithmic quantity we have 

log jdt = 2 log ± log Ani + const, (17) 

where const includes a contribution from the mode mass and 
visibility function. 

The energy supply rate Is believed not to vary with the 
activity level, and the SolarFLAG data were designed to re¬ 
flect this. In other words, we expect 2dA/dX + dF/dX = 0. 
As the d log Ani/dX and d log r„i /dX values are noisy, we 
check for this by performing a weighted least-squares fit with 
the Individual S log Ani and S log Fni as x- and y-values. For 
the 1 = 0,1 values between 2.0 and 3.5 mHz we obtain 
slopes of —1.8 ± 0.58 and —1.47 ± 0.90 for the MLE fitting 
with 100 per centand realistic duty cycle respectively, and 
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-1.53±0.48 and -1.34±0.58 for MCMC. Three of these val¬ 
ues are within la of the expected value of -2, and the other is 
only just outside the Icr range. A fit with d log Ani / dX as the 
abscissa and dlogVni/dX as the ordinate gives agreement 
within the (large) error in all cases, and the mean value of 
2d\ogAni/dX + dlogVni/dX is consistent with zero in all 
four cases. We can therefore say that the fitting is not sig¬ 
nificantly biasing the variation of the energy supply rate, in 
spite of the small bias on the d log Ani/dX and d log Tni / dX 
values. 


4-1.5 Asymmetry 

The design of the SolarFLAG data is such that the fractional 
shift in the asymmetry is expected to be the same as that in 
the line width. Because the asymmetry determinations are 
noisier, we do not show the results per mode. Instead, we 
look at the error-weighted average fractional shift relative 
to the (also error-weighted) temporal mean value for each 
mode. Figure 9 shows the temporal variation of the mean 
fractional asymmetry shift, averaged over all modes between 
2 and 3.5 mHz, compared with that which would have been 
obtained if the fitting returned exactly the input values with 
the same uncertainties. The correlation coefficients between 
the expected and measured mean asymmetry shifts are as 
follows: 0.79 and 0.65 for MLE with 100 per cent and real¬ 
istic duty cycle, and 0.78 and 0.75 for MCMC. The slopes 
of linear least-squares fits with the input average fractional 
shift as the independent variable and the measured values 
as the dependent one are 1.15 ± 0.33, 2.35 ± 0.68 for MLE 
with 100 per cent and realistic duty cycle, and 0.82 ± 0.25, 
1.0 ± 0.49 for MCMC. These results suggest that we can 
(marginally) detect a shift in the asymmetry of about the 
right magnitude and sign, but even when averaging over 
many mode pairs the random errors are too large to make a 
statement about any systematic over- or under-estimate of 
its size, particularly in the realistic-duty-cycle case. We can 
also infer, with caution, that the MCMC fitting is perform¬ 
ing somewhat better than the MLE in this respect. This 
is also true in the case (not shown) where the tests were 
repeated using MLE with the window function handled by 
convolution rather than sidelobe fitting, which suggests that 
the advantage is specifically in the MCMC approach rather 
than the window function convolution. 


4.2 Discussion of the test results 

The results of the simulated-data exercise give us confi¬ 
dence that with 11 years of observations we can determine 
the activity-related shifts in the frequencies without signifi¬ 
cant bias with either fitting method. For the amplitudes and 
line widths both methods may underestimate the shifts by 
around 10 per cent (or more for the 1 = 2 widths). This 
may be the result of a peak profile model that does not fully 
capture the subtleties of the mode asymmetry. We also have 
marginal sensitivity to changes in the asymmetry, even with 
realistic duty cycle. As we are fitting data with asymmetric 
peaks using an asymmetric peak profile, we would not ex¬ 
pect in this case to see any bias in the frequency shifts due 
to neglect of the asymmetry, and indeed, even though we 
only make a weak detection of the asymmetry variation, we 


are very successful in recovering the input frequency shift - 
and also in recovering a zero input frequency shift. A brief 
consideration of the numbers makes this understandable. A 
frequency shift due to neglect of the asymmetry could not 
exceed the product of the line width and the absolute asym¬ 
metry parameter. At 3 mHz the line width is about 1 fiHz 
and the maximum low-degree frequency variation is around 
0.31 /rHz or one-third of the line width, which would require 
around a 30 per cent change in the absolute value of the 
line width - far in excess of what is observed, and so large 
that the asymmetry would be obvious to the eye rather than 
requiring subtle analysis. 

Civen the good agreement between the two techniques 
it seems unlikely that the use of MCMC fitting will inval¬ 
idate the major results from the previous two decades or 
more of MLE fitting of Sun-as-a-star data. 

A comparison of the error estimates returned by the two 
algorithms shows that the uncertainties for the frequency 
are very similar, as are those for the asymmetry parameter. 
For amplitude, we find that the uncertainties from MLE are 
slightly (about 10 per cent) higher than those for MCMC for 
the Z = 0,1 modes and about 15 per cent smaller for Z = 2, 3 
where the signal-to-noise is lower and the modes are made 
up of more components that are handled differently by the 
two algorithms, with the MLE fit having fewer independent 
parameters. For linewidth, the 1 = 0 and I = 1 uncertainties 
from MLE are, not surprisingly, somewhat (20 to 50 per 
cent) larger than those for the 1 = 0/2 and Z = 1/3 pairs 
provided by the MCMC algorithm. 

These findings will guide our interpretation of the re¬ 
sults of the fitting of observational data. 

5 BISON DATA 

5.1 Data 

Early observations by what was to become the BiSON group 
have been made since 1976, but until the early 1990s the 
observations were relatively sparse. The final BiSON net¬ 
work was deployed in 1990-1992 and has been operating 
ever since, although there have been some changes to the 
instruments over time. Fig. 10 shows the yearly duty cycle 
of the network over time. For the purposes of this analysis, 
we use 22 non-overlapping spectra corresponding to the cal¬ 
endar years 1993-2014. We are fortunate in that over this 
period the BiSON duty cycle is essentially uncorrelated with 
the activity level. The data for 1992 and earlier years have 
much lower duty cycle and were excluded from this analysis. 

5.2 Analysis 

5. 2.1 Data preparation 

The BiSON data were prepared as described in Davies et al. 
(2014b) and divided into one-year, non-overlapping time se¬ 
ries. 

5.2.2 Choice of activity proxy 

It is common to express activity changes in terms of a lin¬ 
ear relationship to an activity measure, although sometimes 
additional terms may be used, especially for higher-degree 
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Figure 9. Fitted (symbols) and input (lines) mean fractional asymmetry shift 5b/b as a function of time for MLE (a, c) and MCMC 
(b,d) fits to 11 one-year sets of artificial data, for 100 per cent (a, b) and realistic (c, d) duty cycle. The dashed lines represent the 
activity index scaled to the best linear least-squares fit to the asymmetry-shift data. 


modes. Such analyses have been carried out using, for exam¬ 
ple, the sunspot number, the 10.7 cm radio flux (RF), and 
the global magnetic flux. Bachmann & Brown (1993) found 
that, for observations over the period 1984-1990, the RF and 
the Mgll index gave better correlation with medium-degree 
frequency shifts than the global Kitt Peak magnetic field 
strength, and Chaplin et al. (2007) also found better agree¬ 
ment with proxies such as the RF flux that have a greater 
sensitivity to weak flux. On the other hand, for low-degree 
modes Chaplin et al. (2001) found that the correlations were 
essentially the same for all six proxies they examined, al¬ 
though there was a slight difference between the sensitivity 
to the Kitt Peak magnetic index in the falling phase of Cycle 
22 and the rising phase of Cycle 23. For local and latitudi- 
nally resolved measurements a spatially-resolved proxy is 
needed, and for this purpose a localized measure of the pho- 
tospheric magnetic field strength (often called a magnetic 
activity index, or MAI), tends to be preferred because the 
irradiance-based measures are not spatially resolved. The re¬ 
lationships between these proxies are not necessarily linear, 
and may not be consistent over time as instruments are up¬ 
graded or observing protocols change. Even magnetic mea¬ 
surements from different instruments may show the effects of 
small differences in calibration; also, many magnetic indices 
are only available as Carrington maps that cannot easily be 
translated into daily indices. These considerations need to 
be borne in mind when using the measures to compare with 
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Figure 10. Yearly duty cycle for the BiSON network over the 
period covered by our analysis. 


long sequences of helioseismic data. For the current work, 
as we are dealing with Sun-as-a-Star data, we have chosen 
to use the daily RF index averaged over the exact period 
covered by each spectrum. 
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5.3 Results 

For all variables, we first find the modes that are common 
to all the data sets being compared and then calculate the 
shifts for each mode relative to an error-weighted mean over 
all the data sets. These individual-mode shifts can then be 
combined in an error-weighted mean to show the overall vari¬ 
ation. 

5. 3.1 Frequency 

Figure ll(a,b) shows the mean frequency variation with 
time and with RF index for the two algorithms, and Fig¬ 
ure ll(c,d) shows dUni/dRF for each mode as a function of 
frequency for each method. The methods agree well, even in 
the small deviations from the linear fit; in fact, the corre¬ 
lation coefficient between the MCMC and MLE mean fre¬ 
quency shifts is 0.992, better than the correlation between 
the shifts and the RF index, which is 0.975 in each case. 
The threshold for significance at the 0.1 per cent level is 
0.65 for 22 data points, so all of these are highly signiftcant 
correlations. It is noticeable in Figure 11(a) that in both 
solar minima the frequency points tend to lie below the ac¬ 
tivity trend-line. The frequency dependence of the shifts is 
very clearly seen, and there is a hint that the I = 0 shifts 
are lower than those for I > 0, consistent with the expected 
behaviour due to the latitudinal distribution of the activity 
bands. 

5.3.2 Amplitude 

Figure 12(a,b) shows the amplitude variation for the two 
algorithms as a function of time and as a function of RF 
index, and Figure 12(c,d) shows dlogAni/dRF as a func¬ 
tion of frequency for each algorithm Again, the two methods 
show good agreement, with a correlation coefficient of 0.98 
between the two sets of mean shifts, while the correlation 
with the RF index is -0.88 in each case. 

The data plotted in Fig. 12(a,b) have been corrected 
for the window-function effect on the amplitude, but this 
makes very little difference to the result. The mean shift 
changes by (—0.078 ± 0.005) per cent per RFU from MLE 
and (—0.071 ±0.005) per cent per RFU from MCMC, giving 
a change of about (—8.8T0.6) per cent for MLE or —8.0T0.5 
per cent for MCMC from the highest to the lowest RF value. 
The scatter of the points around the linear 6t, together with 
the agreement between the methods, suggests that effects 
other than a simple linear dependence on activity may be 
involved. Curiously, the points for the years 2011-2013 in 
Figure 12(a) all fall below the ftt to the RF index, suggesting 
that the modes may have been more strongly suppressed in 
the rising phase of Cycle 24. 

The uncertainties are still too large to allow us to use¬ 
fully quantify the frequency dependence of the individual¬ 
mode shifts, but there is a visible tendency for the modes in 
the middle of the five-minute band to be more strongly sup¬ 
pressed by activity. The reduction in obtained by htting 
a quadratic function of frequency to the values with 1 ^2, 
rather than a constant, corresponds to a probability of 12 
per cent for MLE and 23 per cent for MCMC that the result 
could have been obtained from a random distribution - bet¬ 
ter than a la result but not at the 2a level. The shape of the 


variation, with the greatest sensitivity around 3 mHz where 
the modes are strongest, is consistent with the Hndings from 
resolved helioseismology. 

5.3.3 Line width 

Figure 13(a,b) shows the mean shift in logF for the two 
algorithms as a function of time and RF, averaged over 
the modes between 2 and 3.5 mHz. Figure 13(c,d) shows 
dlogF/dRF as a function of frequency for each method. A 
direct comparison between the two data sets in this case is 
more difficult because the MCMC algorithm returns only 
one width per mode pair while the MLE gives two; also the 
different treatment of the window function affects the line 
width estimates. However, the two sets of mean shifts have 
a correlation coefficient of 0.946, while the correlation co¬ 
efficient with the RF index is 0.922 for MLE and 0.951 for 
MCMC . Although the width shifts are larger than those in 
amplitude they also have larger uncertainties. Given the un¬ 
certainties, we do not see any significant deviation from the 
linear trend of the mean shifts with RF index. The mean 
shift is (0.146 ± 0.016) per cent per RFU for MLE and 
(0.142 ±0.013) per cent per RFU for MCMC, corresponding 
to a change of (16.3 ± 1.8) per cent and (15.9 ± 1.5) per cent 
respectively between minimum and maximum RF values. 

With the available data, the statistical significance of 
any frequency dependence of d log Tni/dRF is not high. The 
decrease in obtained by fitting a quadratic function of 
frequency to the values, rather than a constant value, cor¬ 
responds to a 14 per cent probability for MLE and 10 per 
cent for MCMC of this result being obtained from random 
data - again better than a Icr result but not at the 2a level. 

5.3.4 Energy supply rate 

We plot in Figure 14 the variation with date and RF index 
of the mean change in the natural log of the energy sup¬ 
ply rate E„i, formed by taking the sum of the mean line 
width shift and twice the mean amplitude shift. Following 
Chaplin et al. (2000), the « 90 per cent correlation between 
the A and F errors has been taken into account in the error 
propagation. The results imply a fractional change in the 
energy supply rate between minimum and maximum RF of 
(—1.1 ± 2.2) per cent from MLE and (0.17 ± 1.9) per cent 
from MCMC. The previous result that there is no significant 
change in the energy supply rate with activity level is thus 
confirmed. On the other hand, there is a hint of a decrease 
in the energy supply rate from 2011-2013, which we see as 
an unexpectedly large decrease in the mode amplitude in 
this period. 

5.3.5 Asymmetry 

Figure 15 shows the mean fractional shift in the asymmetry 
parameter b as a function of time and of RF index, for the 
two algorithms. As the uncertainties are large, we do not 
show the shifts for individual mode pairs. The coefficient of 
correlation between mean shift and the RF index is 0.36 for 
MLE (just short of being signiftcant at the 10 per cent level) 
and 0.66 for MCMC (significant at the 0.5 per cent level), 
and the correlation between the two sets of mean shifts is 
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Figure 11. Solar-cycle frequency changes from BISON data 1993—2014. The top row shows frequency variation from MLE (filled symbols) 
and MCMC (open symbols) fits, as an error-weighted mean over all modes with I ^ 2 and frequencies between 2 and 3.5 mHz, as a 
function of time (a) and RF index (b). The solid and dashed lines represent linear fits to the RF index for MLE (solid) and MCMC 
(dashed). The bottom row shows du^i/dRF as a function of frequency for MLE (c) and MCMC (d) fits, for I = 0 (black circles), Z = 1 
(blue diamonds), 1 = 2 (red squares) and I = 3 (green open triangles). 


0.46. The results of the artificial data tests lead us to place 
more reliance on the MCMC estimate. Even though the re¬ 
sult is marginally significant, it is consistent with the results 
in Jimenez-Reyes et al. (2007), which only included data on 
the declining phase of the last cycle, where the asymmetry 
change is most pronounced. 


6 DISCUSSION AND CONCLUSIONS 

We have analysed 22 years of BiSON data for activity- 
related changes in the mode parameters using two different 
algorithms. The analysis was validated by the use of an im¬ 
proved process that allowed for the incorporation of the ef¬ 
fects of line asymmetry. The newer MCMC algorithm repre¬ 
sents a small improvement over the older MLE, particularly 
in the search for changes in the asymmetry, but in general 
there is good agreement between the two algorithms. Both 
methods accurately reproduce the input frequency shifts in 
artificial data and slightly underestimate the shifts in mode 
amplitude and width, which is perhaps not surprising given 
the limitations of the model peak profile. 

One of the objectives of the present work was to check 
that the previous reports of frequency shifts were not an 
artefact of underlying changes in the mode asymmetry, as 


suggested by Korzennik (2013). Contrary to this suggestion, 
we find that our fitting accurately reproduces the input fre¬ 
quency shifts in artificial data and finds shifts in the real 
data of the same magnitude as those reported in earlier 
work. Furthermore, the MCMC parameter estimation ex¬ 
plicitly favours a frequency change over a simple change in 
underlying asymmetry (with the model we apply) and in 
any case the asymmetry change required to account for the 
frequency shift would be unrealistically large. 

We have not attempted in this work to reproduce the 
kind of analysis of frequency shifts used by for example 
Basu et al. (2012), in which spectra from heavily overlapped 
time periods were used and the shifts were averaged in fre¬ 
quency bands; our objective was instead to look for system¬ 
atic biases in long-term trends. Having eliminated the fear 
of such biases in the frequency shifts, we will be able to pro¬ 
ceed with confidence to a more detailed analysis of the short¬ 
term variations in further work. We do see deviations from 
the linear relationship between frequency and activity that 
are more evident at low activity and appear consistent with 
those reported by Fletcher et al. (2010); Basu et al. (2012) 
and others. 

Our amplitude and width changes also agree well with 
earlier work. Elsworth et al. (1993) (who htted using a sym¬ 
metric Lorentzian peak profile) report a change of -35 per 
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Figure 12. Solar-cycle mode amplitude changes from fits to 1-year BiSON spectra from 1993-2013. The top row shows amplitude 
variation from MLE (filled symbols) and MCMC (open symbols) fits, averaged over all modes with 1^2 and frequencies between 2 and 
3.5 mHz, as a function of time (a) and RF index (b). The solid and dashed lines represent linear fits to the R,F index for MLF (solid) and 
MCMC (dashed). The bottom row shows shows dlog A^^i/dB as a function of frequency for MLF (c) and MCMC (d), for 1 = 0 (black 
circles), 1 = 1 (blue diamonds), 1 = 2 (red squares) and I = 3 (green open triangles). The dashed curve represents a quadratic function 
in frequency fitted to the modes with 1^2. 


cent in mode area (height times width) between the 1986 
minimum and the 1990 maximum, but based on their linear 
fit to sunspot number a change of -20 per cent would give 
a more realistic comparison with our result. As our ampli¬ 
tude is proportional to the square root of their mode area, 
that would translate to a roughly -10 per cent change in am¬ 
plitude, which compares reasonably well with our estimate 
of around -8 per cent. Chaplin et al. (2000), in an analysis 
covering the period 1991--1997 and averaging over modes be¬ 
tween 2.6 and 3.6 mHz, report changes of (24 ± 3) per cent 
in line width and —(22±3) per cent in power between maxi¬ 
mum and minimum, (corresponding to (—11 ±2) per cent in 
amplitude). The 100-day periods analysed by Chaplin et al. 
(2000) correspond to a range in RF index of about 150 RFU, 
compared with the range of 110 in the 365-day average val¬ 
ues in our analysis, which would make their 24 per cent 
change in line width equivalent to 18 per cent, within a la 
difference from our 15.5 per cent value. For amplitude, their 
change of -11 per cent would be equivalent to (8.1 ± 1.1) per 
cent, in good agreement with our analysis. It is not clear 
whether the « 10 per cent underestimate uncovered by our 
simulations applies to both results; the biases for symmetric 
fits may well be different. 

We do not reproduce the finding of Salabert et al. 


(2007) that the I = 0 mode width is less sensitive to global 
activity than the higher-degree modes. It is not possible to 
check this using the MCMC fitting, as the I = 0 and 1 = 2 
mode components are assumed to have the same width. 

Our results on the asymmetry shifts are somewhat 
disappointing. However, we should note that the work of 
Jimenez-Reyes et al. (2007) was more narrowly focused on 
the asymmetry measurement and they carried out a more 
elaborate analysis to obtain the best possible results. Also, if 
we look only at the eight years of BiSON data that they used 
in the frequency range that they considered, we do reproduce 
their ‘marginally’ significant result for the BiSON data. The 
much weaker correlation that we find for the whole dataset 
suggests that their result for BiSON may have been a statis¬ 
tical fluke; measuring the asymmetry in ground-based data 
with a duty cycle in the 70 per cent range is still challenging. 
In particular, the intimate relationship between background 
and asymmetry may make it difficult to improve very much 
on these results using our current approach of pairwise fit¬ 
ting with a flat background for each pair. We note that the 
asymmetry shifts - and the absolute asymmetry values - 
observed for low-degree modes are too small to have a sig¬ 
nificant influence on the frequency shifts. 

In general, both MLE and MCMC methods perform 
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Figure 13. Width variation with activity. Top row shows mean shift in logF^; from MLE (filled symbols) and MCMC (open symbols) 
fits, averaged over all modes with I ^ 2 and frequencies between 2 and 3.5 mHz, as a function of time (a) and RF index (b). The solid 
and dashed lines represent linear fits to the RF index for MLE (solid) and MCMC (dashed). The bottom row shows dlogV^i/dRF as a 
function of frequency for MLE (c) and MCMC (d) for I = 0 (black circles), I = 1 (blue diamonds), 1 = 2 (red squares), and 1 = 3 (green 
open triangles). For MCMC fits the I = 0 value is for the 1 = 0/2 pair and the I = 1 for the 1 = 1/3 pair. The dashed curve shows a 
quadratic fit to the I ^ 2 values with frequency as the independent variable. 


a b 



Figure 14. Mean change in log energy supply rate from MCMC (open symbols) and MLE (filled symbols) fits, averaged over all modes 
with 1^2 and frequencies between 2 and 3.5 mHz, as a function of time (left) and RF index (right). 


well and give consistent results for the solar-cycle varia¬ 
tions of the mode parameters, as well as very similar uncer¬ 
tainties. This has implications for future peakfinding strate¬ 
gies in both solar and stellar applications. For high-quality 
data the considerable extra computational expense of the 


Bayesian approach may not be justified unless we are specifi¬ 
cally looking for asymmetry changes. We have not addressed 
the issue of the advantages that might be gained by using 
global or pseudo-global instead of pairwise fitting for solar 
data. Such methods are, again, more computationally inten- 
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a 



Date (years) 


b 



Figure 15. Asymmetry variation (mean fractional shift) from MLE (filled symbols) and MCMC (open symbols) fits, averaged over all 
modes with I ^ 2 and frequencies between 2 and 3.5 mHz, as a function of time (left) and RF index (right). The solid and dashed lines 
represent linear fits to the RF index for MLF (dashed)and MCMC (solid). 


sive than the conventional approach, and when combined 
with the Bayesian method this expense becomes nearly pro¬ 
hibitive for large data sets. Our results suggest that it might 
be worthwhile instead to concentrate on developing MLE- 
based global or pseudo-global methods for use with solar 
data; such methods are already in use for asteroseismology of 
solar-like oscillators but may require refinement to deal with 
longer and higher-quality solar data sets. Any such meth¬ 
ods would need to be cross-checked against earlier results. 
Our findings confirm that solar-cycle variations of mode fre¬ 
quency, lifetime, and amplitude from earlier Sun-as-a-star 
work are reliable. If, as in the case of Korzennik (2013), the 
results obtained by novel methods differ substantially from 
those obtained by more conventional means, they should be 
treated with caution. 

The data we have analysed cover the period up to the 
end of 2014, well into the weak but drawn out and double- 
or possibly multiple-humped maximum of Solar Cycle 24; 
we note that the RF value for 2014 is higher than that for 
2013. Some of the issues around inter-cycle differences may 
become clearer once Cycle 24 definitely enters its declining 
phase. We await the next few years of measurements with 
interest. 
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APPENDIX A: THE SOLARFLAG COMPLEX FREQUENCY AMPLITUDE, AND FREQUENCY 
POWER, SPECTRUM 


In the following description - which is based on the detailed discussions in Toutain et al. (2006) Chaplin et al. (2008) - we 
model the p-modes as forced, damped oscillators having a high Q. The frequency response of a given mode (with n, I and m) 
is then just a Lorentzian, which may be written in complex amplitude form as: 

= ^[AnlmO] + (iz)], (Al) 

where: 

' ^nlm 

and 

^1/2 

Q[A^lmC)] = TT^ (^ 2 ) 

' ^nlm 

are the real and imaginary parts, respectively. In the above, Hnim is the height of the peak, and 

^nlm ~ C r^n/m)/(r7iim/2) , 

where Vnim and Fnirn are the central frequency and line width of the mode, respectively. 

The observed response in the frequency domain also depends on the spectral response of the excitation function, which 
we call EnimO- The complex amplitude, AnimO^ nf ^he full response is therefore actually: 

Anlm(i^) ~ LjilmC^Entmiyh (■^^) 


We excite the modes with the granulation-like noise, which has a response that is white only locally in the vicinity of the 
resonance. The power spectral density of this granulation-like noise follows the Harvey (1985) power-law model given in 
Equation 2. To obtain the excitation response, EnimR), we must take the square root of the power spectral density of the 
granulation-like noise; and we must also normalize to the value of the granulation-like noise at the frequency of the resonance. 
This gives: 


Enlm C ) 


/ 1 -I- {2niynlmr)‘^ '\ 
Y 1 + {2 'kut)'^ J 


(A4) 


The power spectral density of this mode is then: 


(A5) 


Now let us write down equations to describe the complete solarFLAG spectrum. This spectrum is comprised of many 
overtones, whose excitation is correlated; and also correlated and uncorrelated background noise. We begin with a description 
of the power spectral density of the overtones of a given (I, m). The excitation of the overtones is correlated, with the coefficient 
p describing this correlation. The modes are also correlated with the granulation-noise background, and this correlation is 
also set equal to p. 

First, consider the equations which describe the real and imaginary parts of the complex amplitude given by the overtones. 
We must sum over all the radial orders n of the chosen {l,m). We therefore have: 






1/2 


over n 


V 


1 -L f2 

' ^nlm 




(A6) 


and 


ilAlmO] 


-E 

over n 


^nlm^nlm 

1 


1/2 


(A7) 


Note that to keep things tidy, we have dropped the explicit dependence of and ^nim on v in the equations above, and in 
what follows. The observed power spectral density of the overtones and the granulation-like noise for this (Z, m) is then given 
by: 


where 


PSDi„i(j/) = IpI tplmOi’lmO + (1 - IpI) 



Hr,l m Enl m 

1 + Cnim 



= ^[AlmO] -b i3[Ai^(i/)] - . 


(A8) 


(A9) 


There are two important things to notice about Equations A8 and A9. First the description in Equation A8 has two parts: 
one part, prefixed by the factor \p\, describes the correlated contributions to the spectrum; while the other part, prefixed by 
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the factor (1 — |p|), describes the uncorrelated contributions. Second, notice how the background factor in Equation A9 
is prefixed by a minus sign. This makes the correlation of the modes and background negative, which in turn means we get 
negative asymmetry as in real Doppler velocity observations. 

Finally, the total power spectrum density is given by summing, incoherently, the power spectral densities of all {I, m) in 
the spectrum, and the uncorrelated photon shot-noise background, Upsn (see Section 2.2.5), to give: 

PSD(i/) = ^ PSDi,„(i/) + Upsn. (AlO) 

all {l,m) 
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